Mott transition in bosonic systems: Insights from the variational approach 
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We study the Mott transition occurring for bosonic Hubbard models in one, two, and three spatial 
dimensions, by means of a variational wave function benchmarked by Green's function Monte Carlo 
calculations. We show that a very accurate variational wave function, constructed by applying a 
long-range Jastrow factor to the non-interacting boson ground state, can describe the superfluid- 
insulator transition in any dimensionality. Moreover, by mapping the quantum averages over such 
a wave function into the the partition function of a classical model, important insights into the 
insulating phase are uncovered. Finally, the evidence in favor of anomalous scenarios for the Mott 
transition in two dimensions are reported whenever additional long-range repulsive interactions are 
added to the Hamiltonian. 

PACS numbers: 74.20.Mn, 71.10.Fd, 71.10.Pm, 71.27.-|-a 
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I. INTRODUCTION 

The recent advances achieved on cold atoms trapped 
in optical lattices have generated an increasing interest in 
the condensed matter community, since they allow exper- 
imental realizations of simple lattice models. ^'^ A great 
advantage of these systems is the possibility to have a 
direct control of the parameters, such as the width of the 
bands and the strength of the interactions, that can be 
manipulated by varying the depth of the optical poten- 
tial.'^ Therefore, cold atoms on optical lattices give the 
unique opportunity to make a close contact with theoret- 
ical models and to examine the origin of the fundamen- 
tal physical phenomena that occur in crystalline materi- 
als. In particular, one of the most spectacular example is 
given by the superfluid-insulator transition in a system 
of interacting bosons: The so-called Mott transition.^ 

In this paper, we consider the superfluid-insulator 
transition in a bosonic system with different kinds of in- 
teractions. The simplest model that contains the basic 
ingredients of strong correlations is the Hubbard model 



h.c. 



1), (1) 



where (...) indicates nearest-neighbor sites, 6| (bi) cre- 
ates (destroys) a boson on the site i, and rii = blbi is the 
local density operator. Here, we consider N particles on 
a one-dimensional (ID) chain, a two-dimensional (2D) 
square lattice, and a three-dimensional (3D) cubic lat- 
tice with L sites and periodic-boundary conditions. At 
zero temperature and integer densities p — N/L there 
is a superfluid-insulator transition when the ratio be- 
tween the kinetic energy and the on-site interaction is 
varied. Otherwise, for non-integer fillings, the ground 
state is always superfluid. In a seminal paper,* by using 
a field-theoretical approach. Fisher and coworkers pro- 
posed that the transition of the d-dimensional clean sys- 
tem belongs to the XY universality class in d + 1. This 



scenario has been confirmed mostly in one and two di- 
mensions by using different numerical techniques, such 
as quantum Monte Carlo and density-matrix renormal- 
ization group. ^i^i^i^i^ In particular, it has been verified 
that in one dimension at p = 1 there is a Kosterlitz- 
Thouless transition and the estimation of the critical 
value of the on-site interaction ranges between Uc/t ^1.7 
and Uc/t ^ 2.3.— Instead, a second-order phase transi- 
tion is claimed to occur for Uc/t ~ 8.5 in two dimen- 
sions?^ and for Uc/t ^ 14.7 in three dimensions J^. Ac- 
curate estimations for the critical values can be also ob- 
tained from a strong-coupling expansion.— 

Besides the numerically exact techniques, important 
insights into the various phases can be obtained by con- 
sidering simplified variational wave functions. The sim- 
plest example is given by the celebrated Gutzwillcr state, 
where the on-site correlation term allows for the suppres- 
sion of the energetically expensive charge fluctuations. 
Contrary to the fermionic case, when considering bosons, 
it is possible to deal with this wave function without any 
further approximation.-i^iM, Then, the Mott transition is 
obtained with a reasonable estimate of the corresponding 
critical value Uc/t, namely Uc/t — d{^/ri^-\-y/nc -\- 1)^ for 
integer fillings p = Uc- The main drawback of this ap- 
proach is that, similarly to what happens with fermions, 
the transition is reached with a vanishing kinetic energy 
and the insulating state completely lacks charge fiuctua- 
tions, namely all particles are localized, Uc in each lattice 
site. Of course, this leads to a wrong description of the 
insulator, whenever the local interaction is finite. 

Following the ideas of previous works on fermionic sys- 
temsfi^ii^iii in a recent paper^* we have shown that, in 
order to correct this outcome and obtain a more accu- 
rate description of the ground state, it is necessary to 
include a long-range Jastrow factor, whose singular be- 
havior at small momenta was shown to be able to turn a 
non-interacting bosonic state into an insulator that still 
contains density fluctuations. In this paper, we present 
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a more detailed study of the properties of that Jastrow 
wave function, as well as of its accuracy in comparison 
with Green's function Monte Carlo (GFMC) calculations, 
which, because of the absence of the sign problem, pro- 
vide numerically exact resultsJ-?'^*^ 

The same Jastrow wave function turns out to be also 
very effective to describe Hamiltonians that contain long- 
range interactions. This case has not been considered 
much in the literature, and we will show that different 
scenarios for the superfluid-insulator transition could oc- 
cur. As a matter of fact, long-range interactions have 
been studied mostly in the continuum, where a transition 
between a charged bosonic fluid and a Wigner crystal has 
been found by varying the density— i^Sii^^ In particular, in 
the presence of a logarithmic repulsion, the 2D Bose liq- 
uid was found to have no condensate fraction, due to the 
predominance of long- wavelength plasmon excitations.— 
In the last part of this paper, we generalize the Hubbard 
model of Eq. (P) to 



condensed into the q = state 



(2) 



where il,{Ri,Rj) is a long-range potential that only de- 
pends upon the relative distance \Ri — Rj\ between two 
particles and V represents its strength. In particular, 
we will consider two possibilities for the long-range po- 
tential. The first one is obtained by taking the Coulomb 
interaction between (charged) bosons moving in a 2D lat- 
tice embedded in a three-dimensional environment, which 
leads to the following potential in reciprocal space 



1 



(3) 



y/ {cosqx + cosqy - 3)^ 

with a small-q behavior fl{q) ^ l/l?!- The second case 
consists in directly considering the solution of the Poisson 
equation in 2D: 

^iqx,qy) = 7, -, ^ W 

2 — [cosqx + cosqy) 

which for small momenta behaves like ri(q) ~ l/?^, 
leading to a logarithmic interaction in real space, i.e., 
r2(r) ~ — log(r). In both uniform background is 

considered in order to cancel the divergent q = Q compo- 
nent of the potential. 

The paper is organized as follow: In section [Til we de- 
scribe the variational wave function, in sections lllll IIVI 
andfVl we show the results for the ID, 2D and 3D short- 
range models. Then, in section IVIl we consider the 
2D case with long-range interactions and finally, in sec- 
tion |VlIl we draw our conclusions. 



II. THE VARIATIONAL WAVE FUNCTION 
A. The Jastrow wave function 

The variational wave function is defined by applying 
a density Jastrow factor to a state with all the bosons 



exp ■ 



2 ^ 



Vi,j{ni ~ l)inj 



\^o), (5) 



where |$o) ~ the non-interacting boson 

ground state with N particles, (n^ — 1) is the variation 
of the on-site density with respect to the average value 
p = 1, and Vij are translationally invariant parameters 
that can be optimized to minimize the variational en- 
ergy*^ In order to get some physical insight from the 
variational state, it is more instructive to consider the 
Fourier transform of the Jastrow parameters Vq. Indeed, 
there is a tight connection between the small-q behavior 
of Vq and the nature of the ground state. In particular, 
as we are going to discuss, Vq ^ ^/\q\ implies the exis- 
tence of sound modes in any dimensions, as expected in a 
superfluid. On the contrary, to recover an insulating be- 
havior a much more singular Vq for g — s- is required 
The physical reason is that, in order to describe a re- 
alistic Mott insulating wave function that does include 
charge fluctuations, it is necessary to spatially correlate 
the latters. This is accomplished by a sufficiently sin- 
gular Vq that favors configurations where opposite-sign 
fluctuations, (n^ — l){nj — 1) < 0, are close to each other, 
while equal-sign ones, (n^ — l)(nj — 1) > 0, are far apart. 

It should be mentioned that previous studies of 
fermionic systems^ and more recent ones on the bosonic 
Hubbard modelf^^ stressed the importance of a many- 
body term containing holon-doublon interactions for 
nearest-neighbor sites: 

|*g,A/s> =exp (^~gJ2^1 +9MBj2i}j |*o), (6) 

where g and gMB are variational parameters and the 
many-body operator is defined by 

i^ = Y[{1 - d,+s) + ck Hil - h,+s), (7) 



where hi — 1 (d^ = 1) if the site i is empty (doubly occu- 
pied) and otherwise, S = zLx,zLy; therefore, counts 
the number of isolated holons (empty sites) and doublons 
(doubly occupied sites). Even though this operator has 
been originally introduced for fermionic systems, where 
the maximum occupancy at each site is given by two 
electrons, it is useful also for bosons, since in the limit of 
large U/t the number of sites with an occupation larger 
than two is negligible. However, within this framework, 
the evidence for a true Mott transition is rather contro- 
versial, and it is not clear if an insulating phase can be 
stabilized in the thermodynamic limit. 

Combining together the variational Eqs. ^ and 
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we obtain the variational ansatz 



1^ 



jmb) = cxp 



{n, - l){nj - 1) 



(8) 



containing both the long-range Jastrow factor and a 
short-range many-body term. As it will be shown in the 
next sections, the presence of the latter term is impor- 
tant in 2D and 3D, mainly to increase the accuracy in the 
strong-coupling regime. Instead, in ID, the many-body 
term does not improve the accuracy of the long-range 
Jastrow state and there is no appreciable difference be- 
tween the wave function ([5]) and ([5]). In any case, we 
emphasize that, according to our calculations, the short- 
range term alone cannot lead to an insulating behavior, 
and the main ingredient to drive a superfluid-insulator 
transition is the long-range Jastrow factor, parametrized 

by Vq. 



B. Mapping onto a classical model 

The variational calculation with the wave function ([S]) 
can be shown to correspond to a classical problem at fi- 
nite temperature. This correspondence provides many 
insights into the properties of \'^j,mb)- To prove the 
mapping, let us denote a bosonic configuration by the po- 
sitions {x} of the particles and a generic quantum state 
by 1^'). For all the operators 9 diagonal in space coordi- 
nates, the quantum average 



(9) 



can be written in terms of the classical distribution 



|vi/(x)p = |(x|f)|VE.'l(^'l*)P,as 

{9)=J2{x\9\x)\^ixW 



(10) 



Since |5'(a;)p is a positive quantity, there is a precise cor- 
respondence between the wave function and an effective 
classical potential Vci{x): 



\^ix)\ 



2 ^ e-Ki(x)_ 



(11) 



The explicit form of the potential Vci{x) depends upon 
the choice of the Jastrow factor, while |$o) does not 
contribute to it, since (a;|$o) does not depend upon the 
configuration \x). In particular, whenever there is only 
the two-body potential (i.e., gMB — 0) Vc\{x) is Gaus- 
sian, i.e., Vciix) = J2q=io^'i'^<!(^)'^-(i(^)^ being nq{x) the 
Fourier transform of the local density of the configuration 
\x). In this case, the variational problem becomes equiva- 
lent to solve a classical model of oppositely charged parti- 
cles (the holons and the doublons) mutually interacting 



through a potential determined by Vq. In the presence 
of the short-range many-body term, V^i(x) is no longer 
Gaussian. However, when the density fluctuations are 
suppressed at large U/t, the quadratic term gives the 
most relevant contribution, hence the mapping onto a 
classical model of interacting oppositely charged parti- 
cle still holds with Vc\{x) ~ Tlq^a'^'c/ ^ ''^q{^)''^-q{^)^ 
being replaced by a slightly different effective potential 
Vq^^ . In spite of the differences, it is plausible that the 
small-(j behavior of Vq-^-^ must follow the same singular 
behavior of the Jastrow potential Vq. 

Let us now discuss in more details the connection be- 
tween the form of Vq and the low-energy excitation spec- 
trum of the system. By means of the /-sum rule one can 
show thalii^iSl 



, , D I du;ujN{q,uj) 

i?„ = -24^^sin^ (|) = 



, (12) 



where N{q,Lu) is the dynamical structure factor, Nq = 
/ duj N{q, to) the static one, and (k) the hopping energy 
per site. Eq can be interpreted as the average excitation 
energy of density fluctuations at momentum q. Eq 
for q —> is a sufHcient but not necessary condition for 
the existence of gapless excitations. In particular, deep 
inside the superfluid phase, Eq must coincide with the 
energy dispersion of the Bogoliubov sound mode, which 
carries most of the spectral weight. Although we can not 
access directly dynamical properties by our variational 
wave function, still we can provide a variational estimate 
of Eq, in the same spirit of the Feynmann's construction 
in liquid Helium This amounts to use the variational 
values of (k) and of the static structure factor Nq, defined 
as 



(13) 



with Uq the Fourier transform of the local density. Note 
that the uncorrelated 



($o|n_^nq|^o) 
(•Sol^o) 



(14) 



is constant for any q ^ at p = 1. Whenever a weak- 
coupling approach in the Jastrow potential Vq is possible, 
the following relation holds 



q 



1 



2NOVq 



Ar^ = ^ng(a:)n„,(a:)e-^='(") 

(15) 

the last equality following from the singular behavior of 
Vq that is expected both in the superfluid and in the 
Mott insulating phases. Eq. (|15p shows that Vq ^ l/l^l 
allows to recover the correct behavior oi Nq ~ \q\ and 
Eq \q\ in the superfluid regime, which is what we in- 
deed find, see following sections. In the insulating phase, 
should Eq. (fT5|) be valid, we would expect Vq ~ P/q^ 
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to get the expected behavior Nq ^ and Eq finite for 
small q. This would correspond through Eq. pT|) to a 
classical Coulomb gas model (CGM) with effective tem- 
perature Teff — IT / (3^- In ID, for any finite temperature, 
the CGM is always in a confined phase where oppositely 
charged particles are tightly bound in pairs, and with 
exponential decaying correlation functionsi^ This sug- 
gests that the ID CGM may indeed provide through the 
mapping pT|) a good description of a ID Mott insulating 
wave function. 

In 2D the CGM undergoes a Berezinskii-Kosterlitz- 
Thouless phase transition at T^'^^ ^ 1/4, between a 
confined phase (stable at low temperature) and a plasma 
phase (stable at high temperature). Similarly to the ID 
case, one would argue that the confined phase should 
correspond to the 2D Mott insulator. However, the 2D 
confined phase of the CGM displays power-law decaying 
correlations'^*^ that would correspond to power-law decay- 
ing equal-time correlations of the quantum ground state. 
This is not compatible with a genuine Mott insulator, 
which already suggests that the insulating wave function 
must be characterized by a Jastrow potential Vq more 
singular than 1/q^ as g — > 0, as indeed we find. In turns, 
this implies that (jlSp must not be valid, since, in spite of 
Vq X q'^ ^ oo for q ^ 0, we still expect Nq ^ q^ . 

The breakdown of Eq. (jlSp becomes even more pro- 
nounced in 3D, where a potential Vq ^ cannot de- 
scribe at all an insulator since it is not sufficiently sin- 
gular to empty the condensate fraction. ■^^ In fact, we 
find that the optimized variational wave function shows a 
more diverging Vq ^ in the 3D insulator though 

Nq ~ q^ . The properties connected to the insulating 
phase can be again uncovered within a classical 3D gas 
with a l/I^P potential, recently considered in Ref. [SJ. 
Indeed, analogously to what happens in 2D, also in 3D a 
system of charges interacting with a logarithmic poten- 
tial in real space admits a high-temperature fiuid regime 
separated by a low-temperature dielectric phase. Within 
this mapping, the 3D insulating state is found to cor- 
respond to the low-temperature phase of this classical 
model. 

In the following, we will present our results obtained by 
considering the quantum variational wave function and 
we will use the classical mapping to gain insights into the 
ground-state properties. Moreover, in order to verify the 
accuracy of the variational calculations, we will perform 
the numerically exact GFMC that allows one to obtain 
ground-state properties. 



III. THE ID HUBBARD MODEL 

Here, we consider a chain of L sites with periodic 
boundary conditions and N = L bosons. First of all, in 
Fig. [T] we compare the variational accuracy of the wave 
functions ^ and ([5]) for different values of U /t. Once 
a long-range Jastrow factor is considered, the many-body 
term of Eq. ([7]) is irrelevant and there is no appreciable 




12 3 4 

U/t 

FIG. 1: Accuracy of different variational wave functions as a 
function oiU/t for 60 sites and 60 bosons. AE = Eq — Evmc , 
where Evmc is the variational energy and Eo is the ground- 
state one, obtained by GFMC. The state of Eq. ([5} is denoted 
by "Jastrow" , the one of Eq. (O by "Gutzwiller-|-MB" , and 
the one of Eq. dS} by "Jastrow-|-MB" . 
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FIG. 2; Jastrow parameters Vq multiplied by as a function 
of \q\ for 60 sites and 60 bosons. 



difference between the wave functions ([5]) and ([8]), for all 
the values of the on-site interaction. By contrast, the 
Gutzwiller state, also when supplied by the many-body 
term, is much less accurate by increasing U jt. 

Therefore, in the following, we will consider the state 
with long-range Jastrow factor ([5]) alone, since the many- 
body term makes the algorithm much slower and does not 
improve the quality of the variational state. In Fig [21 we 
report the minimized Jastrow parameters Vq multiplied 
by <f for different U jt: There is a clear difference in the 
small-q behavior for U jt < 2.4, where Vq ~ ct/\q\ and for 
U/t > 2.5, where Vq ~ /3/q'^. At the variational level, the 
change of the singular behavior of the Jastrow parame- 
ters for U/t ~ 2.45 marks the superfluid-insulator tran- 
sition. Indeed, as discussed in the previous paragraph, 
Vq ~ ct/\q\ implies a gapless system, whereas Vq ~ P/q"^ 
indicates a finite gap in the excitation spectrum. Let 
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FIG. 3: Density structure factor A^g divided by |g| calculated 
with variational Monte Carlo (left panel) and GFMC (right 
panel) for different U/t and L = 60. From top to bottom 
U/t = 1.6, 1.8, 2, 2.2, 2.4, 2.5, 3, and 4. 



FIG. 4: Variational results for the momentum distribution n^! 
in ID for L = 60 (squares), 100 (circles), and 150 (triangles) 
across the transition {U/t = 2.4 and 2.5). Inset: Size scaling 
of the condensate fraction no/L for U/t = 2.4 (upper curve) 
and U/t — 2.5 (lower curve). 



us now concentrate on the insulating phase. Here, the 
Jastrow wave function ^ can be mapped onto the par- 
tition function of an effective classical CGM and f3 plays 
the role of the inverse classical temperature /3 = ir/Teff. 
In ID the CGM is in the confined phase for any finite 
temperature, with exponential correlations This out- 
come is consistent with the fact of having, in the quantum 
model, a finite gap in the excitation spectrum. Remark- 
ably close to the Mott transition in the insulating phase, 
the value of (3 obtained from the optimized Jastrow po- 
tential is very small and approaches to zero when U Uc 
from above. Since the correlation length of the ID CGM 
diverges for /? — > 0, our numerical outcome gives a strong 
indication in favor of a continuous transition between the 
superfluid and the Mott insulating phase. 

Let us now analyze the transition by considering the 
density structure factor Nq. In the small-q regime, we 
can generally write that 



^9 =7ikl+72'z' + 0(g3), 



(16) 



where 71 and 72 depend upon the Jastrow parameters. 
In analogy with spin systems, we have that 71 = VcX, 
with Vc and x the sound velocity and the compressibility, 
respectively. The fact of having 71 = in the insulating 
regime indicates that this state is incompressible (i.e., 
X = 0). From Fig. [3l we obtain that 71 has a very sharp 
crossover from a finite value to zero across the transi- 
tion, suggestive of a true jump in the thermodynamic 
limit. This outcome is consistent with the fact that the 
compressibility also has a finite jump in the ID quantum 
phase transitionj2i Moreover, just above Uc in the insu- 
lating regime, 72 is very large for both the variational and 
the GFMC calculations, indicating the peculiar character 
of the ID transition. By using the small-g behavior of the 
density structure factor, the Mott transition can be lo- 
cated at Uc/t ^ 2.45±0.05 for the variational calculation, 
whereas the GFMC approach gives Uc/t ~ 2.1 ± 0.1. 



0.8 
0.6 
0.4 
0.2 







♦ 


— ♦ ♦ — 

t tt 


— ^ 
















♦ 







0.02 

1/L 



2 
U/t 



FIG. 5: Superfluid stiffness Da calculated by GFMC as a 
function of U/t for different sizes in the ID boson Hubbard 
model. Lower inset: Size scaling of Ds for different U/t (same 
values of the main panel, with U/t increasing from top to 
bottom). Upper inset: Ds x L as a function of U/t. The 
point where the different curves cross marks the transition 
point. 



The superfluid-insulator transition can be also easily 
detected by considering the momentum distribution: 



{^j\bW\^j) 



J) 



(17) 



where b\, is the creation operator of a boson of momentum 
k. This quantity has a radically different behavior below 
and above the transition: In the superfluid phase, it has 
a cusp at fc = 0, although there is no condensate fraction, 
i.e., hq/L — > for L 00, while in the insulating phase 
it is a smooth function of the momentum fc, see Fig. 
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Finally, we want to conclude the ID part by consider- 
ing the superfluid stiffness Dg. In analogy to what has 
been done by Pollock and Ceperley at finite tempera- 
ture^^ this quantity can be also calculated directly at 
zero temperature by using the GFMC and the so-called 
winding numbers (see Appendix) 



Ds = lim 



(vI^oIIW^(t)PIvI^o) 
D Lt 



(18) 



where |5'o) is the ground-state wave function obtained by 
GFMC, D is the dimension of the system and W{t) = 
^Jri(T) — ■ri(O)], ri(T) being the position of the i-th par- 
ticle after evolving it for a diffusion time r from the initial 
position ri(0). The diffusion process must be done with- 
out considering periodic boundary conditions, namely by 
increasing or decreasing the values of the coordinates of 
a particle that crosses the boundaries of the cluster. In 
this way, non-zero winding numbers across the lattice 
can be detected. It should be stressed that, exactly at 
zero temperature, Dg can only give information about 
the presence of a gap in the excitation spectrum, and, 
therefore, it can discriminate between conducting and 
insulating phases. Our results show that Dg is finite and 
large in the weak-coupling regime, whereas it vanishes 
for U /t > 2.1, see Fig. [SJ Again, in analogy with spin 
systems and from general scaling arguments valid for ID 
boson models, we expect a jump at the transition, that 
however is very hard to detect by considering finite clus- 
ters;^ Nevertheless, an accurate value of the transition 
point is obtained from the size scaling of Dg x L, see 
Fig. El 



IV. THE 2D HUBBARD MODEL 

Let us now turn to the 2D Hubbard model and consider 
square clusters with L = I x I sites and N = L bosons. 
The accuracy of the three wave functions ([5]) , (O and ^ 
are reported in Fig. [6] The situation is different from the 
previous ID case: The Gutzwiller state with the many- 
body term, which in ID is not accurate for large U/t, in 
2D becomes competitive with the Jastrow wave function 
for U/t> 14. Moreover, the presence of the many-body 
term strongly improves the accuracy of the Jastrow state 
as soon as U/t > 8. Then, in the following, we will 
consider the state ([8]) for both the variational and the 
GFMC calculations. 

In Fig. [71 we show the behavior of the optimized Vq as 
a function of the interaction strength. Similarly to what 
happens in the ID system, we obtain that Vq ~ o:/\q\ for 
U/t < 10.5, while Vq is best fitted by ~ — \og{q)/q^ for 
U/t > 10.5 (see Fig. [5]), corresponding to the superfiuid 
and the Mott insulating phase, respectively. As we an- 
ticipated in section fll Bl the Jastrow potential in the in- 
sulating phase is more singular than 1/q^, suggestive of a 
classical model with bound charges but presumably with- 
out the power-law correlations displayed by the CGM 
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FIG. 6: Accuracy of different variational wave functions as 
function of U/t for the 10 x 10 cluster and 100 bosons. The 
symbols are the same as in Fig. [Tj 




FIG. 7; Jastrow parameters Vq multiplied by as a function 
of |g| [along the (1,0) direction] for 20 x 20 (circles), 26 x 26 
(squares), and 30 x 30 (triangles) clusters. Inset: The many- 
body variational parameter qmb as a function of U/t. 



in the confined phase. In reality, on the sizes available 
to our numerical approach, we can not firmly establish 
whether a classical potential —\og{q)/q'^ has indeed ex- 
ponential decaying correlation functions. Nevertheless, it 
is remarkable and very encouraging that the variational 
optimization leads to such a singular Vq. 

The evidence that the change in behavior of Vq does 
correspond to the superfluid-insulator transition comes 
also by the small-g limit of the structure factor Nq. In 
Fig. [9l we show the results for the variational and the 
GFMC calculations as a function of U/t. In both cases, 
we find a different small-g behavior for large and small 
couplings. In the variational calculations, for U/t < 10.3 
the structure factor behaves as Nq ^ "yi\q\, while for 
U/t > 10.3 we get Nq ~ 729^- Therefore, at the varia- 
tional level, the superfiuid-insulator transition is located 
around Uc/t ^ 10.3 ± 0.1; this value is slightly smaller 
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FIG. 8: Jastrow parameters Vq multiplied by g^/ log((j'/7r) as 
a function of l^l [along the (1,0) direction] for different U/t 
and the same sizes as Fig. [7] 



FIG. 10; GFMC results for the sound velocity Vc obtained 
through a finite size scaling of the ground-state energy, see 
Eq. (|19p . The line is a guide to the eyes. 
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FIG. 9: Left panel: Density structure factor A^, divided by 
|g| calculated by the variational Monte Carlo for different U/t 
and L = 20 X 20. From top to bottom U/t = 10, 10.2, 10.4, 
11, and 12. Right panel: The same for the GFMC on the 
i = 16 X 16 cluster. From top to bottom U/t = 8, 8.2, 8.4, 
8.6, and 8.8. 
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FIG. 11: Variational results for the momentum distribution 
Rfc in 2D for 10 x 10 (squares), 16 x 16 (circles), 20 x 20 
(upward triangles), and 26 x 26 (downward triangles) clusters 
with U/t = 10 and 12. Inset: Size scaling of the condensate 
fraction no/L for U/t = 10 (upper curve) and U/t = 12 (lower 
curve) . 



than the one obtained by the simple Gutzwiller wave 
function, for which Uc/t ~ 11.65.— The critical value of 
the on-site interaction is rather different within GFMC, 
yielding Uc/t ^ 8.5±0.1, in close agreement with the 
value found in the literaturci^ Let us focus more deeply 
on the behavior of the structure factor Nq. Approach- 
ing the transition from the weak-coupling region, 71 goes 
smoothly to zero, in contrast with the results of the ID 
model, where we observed an abrupt jump. Also con- 
trasting the ID case, 72 is found to be finite and continu- 
ous across the transition. These results, based upon the 
variational wave function, are confirmed by the GFMC 
calculations, see Fig. [9l The vanishing linear coefficient 
of Nq can be ascribed either to Vg or to x- In order 
to clarify which one of these quantities goes to zero at 
the transition, we extract the sound velocity Vs from the 



finite-size scaling of the exact ground-state energy by 

6o(i) = eo(oo)-^, (19) 

where e.o{L) is the ground-state energy per site for a clus- 
ter with L = P sites, eo(oo) is the extrapolated value 
in the thermodynamic limit, and cq is a given model- 
dependent constant. Our results, shown in Fig. I10[ 
clearly indicate that stays finite across the superfluid- 
insulator transition, thus implying a vanishing compress- 
ibility when approaching the insulator. 

The fingerprint of this transition is also given by the 
momentum distribution, see Fig. 1111 For this quantity, a 
striking difference is observed below and above Uc- In the 
former case, a cusp- like behavior with a finite condensate 
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a function of U/t for the 10 x 10 x 10 cluster and 1000 bosons. 
The symbols are the same as in Fig. [T] 



fraction is found, whereas in the latter case a smooth 
behavior is detected, with a vanishing ilq/L. Notice that 
a vanishing condensate fraction in the thermodynamic 



hmit immediately follows from Nq 
sum rule derived in Ref. l36i . 



q by using the f- 



Finally, we report in Fig. [T^the stiffness Dg , calculated 
by using GFMC. In 2D, we obtain a different behavior 
with respect to the ID case, where a finite jump is rather 
clear at the superfluid-insulator transition. Indeed, the 
evaluation of the stiffness for different sizes confirms the 
absence of the jump in 2D, as expected for a second order 
phase transition!^ 
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for non-optimized wave functions with Vq ~ /^ao/lgp for two 
values of Pzd- 



V. THE 3D HUBBARD MODEL 

Let us move now to the 3D case. Here, we mostly re- 
strict our analysis to the variational method that allows 
us to assess rather large sizes. We start by considering the 
accuracy of the different variational states as a function 
of the interaction J7/t, see Fig. [131 It turns out that in 
the large-?/ region the most accurate state contains both 
the long-range Jastrow term and the many-body term, 
as it occurs in 2D. Moreover, in analogy to ID and 2D, 
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the presence of a phase transition upon increasing U /t 
is clearly signaled by the sudden change in the small-^/ 
behavior of the Jastrow factor, see Fig. [Ml Here, its 
behavior changes drastically from Vq ^ (x/\q\ to a more 
diverging Vq ~ l/I^P in the insulating regime. In partic- 
ular, the sudden change of behavior allows us to locate 
the transition around Uc/t ~ 18, which is close to the 
critical value of recent Monte Carlo simulations.^" 

Focusing on the large-t/ region of the phase diagram, 
let us discuss the implications of the singular Vq ^ 
Jastrow potential. First of all, following the same argu- 
ments of Ref. HH, we can assert that the strong singular 
character of the Jastrow potential is able to empty the 
condensate, whereas a less singular Vq ^ 1 jq^ would not 
lead to no — > in the thermodynamic limit. Remarkably, 
even though Vq ^ l/jqp, the structure factor in the in- 
sulator has the correct Nq 
turn, this implies that Eq. 
prove more firmly that Vq 



q' behavior, see Fig. [TS] In 



15)) does not hold. In order to 
^ Psd/IqI'^ can indeed lead to 



J3D/ 

Nq^ q^, we have calculated the density structure factor 
with a non- optimized wave function of the form ([5]) with 
Vq ~ /Ssd/I^P and for different values of Psd- As shown 
in Fig. [iSl for small (Ssd we obtain Nq ~ \q\^^, implying 
that Eq. is qualitatively correct in this case. How- 
ever, above a critical P^jj, the small-q behavior of the 
density structure factor turns into Nq ^ g^, signaling a 



remarkable breakdown of Eq. (fT5)) . From Fig. [T3]it turns 
out that the optimal value of Psd that we get variation- 
ally at the superfluid-insulator transition is larger than 
f3^jj, confirming that Nq ~ q^ in the insulating phase. 
Most importantly, the change of behavior as a function 
of PsD is consistent with the binding-unbinding phase- 
transition recently uncovered in a classical 3D gas with a 
l/|gp potential. ^■^ Once again, as in ID and 2D, the Mott 
insulating wave function in 3D is found to be closely re- 
lated to the low-temperature confined phase of a classical 
model, where opposite charges tightly bound. 
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VI. 2D SYSTEMS WITH LONG-RANGE 
INTERACTION 



FIG. 17: Variational results for the density structure factor 
N{q) divided by \qf^'^ for the potential of Eq. Q with differ- 
ent V/t and 20 x 20 (circles), 26 x 26 (squares), and 30 x 30 
(triangles) clusters. 



In the previous paragraphs, we have shown that, in 
the Hubbard model H]), the superfluid regime can be 
described by a long-range Jastrow wave function with 
Vq ~ ct/l*?!- By increasing the on-site interaction U, 
our variational approach describes a continuous transi- 
tion to an insulating phase that corresponds to the con- 
fined phase of a classical model of interacting particles 
with opposite charge. In particular, we have found ev- 
idences that the 2D Mott insulating wave function cor- 
responds to a classical model with —\og{q)/q^ potential 
rather than to a 2D CGM with potential P/q'^ in the 
confined phase {P > (3* > 47r). This result, as we dis- 
cussed, has a physical importance since the confined 2D 
CGM has power-law correlations, not possible in a Mott 
insulating phase. Nevertheless, it would be interesting to 
search for bosonic Hamiltonians whose ground state can 



be described by the variational wave function ([8|) with a 
Jastrow potential Vq ~ P/q^ ■ 

For that purpose, we consider the more general Hub- 
bard model of Eq. ([2]) with a long-range interaction il(r). 
Let us start by considering the realistic case of a Coulomb 
potential, namely r2(r) ~ l/r that in 2D leads to Eq. ([3]). 
Then, we can vary its strength V to drive the system 
across a superfluid-insulator transition. The small-g be- 
havior of the optimized Jastrow parameters is shown in 
Fig. [m For V/t < 8, we obtain that Vq ~ l/|gp/^. 
On the contrary, for larger values of the interactions, Vq 
becomes more singular and in particular is best fitted by 
— \og{q) / q^ , just like what we found for short-range inter- 
action. The structure factor in the weak-coupling phase 
behaves as Nq ~ |(7p^^, see Fig. [T71 compatible with a 
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FIG. 18: The same as in Fig. [HI for the potential of Eq. (gj. 



/3 ^ /3* at the transition. This indicates that the CGM 
that corresponds to the variational wave function is in the 
plasma phase, with exponential decaying density-density 
correlation functions. Indeed, the structure factor be- 
haves like Nq ~ for all the coupling strengths V/t, see 
Fig. [ini We believe that the weak-coupling phase has to 
be identified with the algebraic long-range ordered phase 
found at high density in the continuum model by Magro 
and Ceperleyi^ This phase is characterized by absence 
of condensate but by a power-law decay of the single- 
particle density matrix. On the contrary, the strong- 
coupling phase with Jastrow potential Vq ^ — log(g)/g^ 
must correspond to a genuine Mott insulator with all cor- 
relation functions decaying exponentially. 



VII. CONCLUSIONS 



- V/t=6 




|q| 

FIG. 19: The same as in Fig.[T7l for the potential of Eq. Q. 

superfuid phase with 2D plasmons. We mention that 
a similar behavior has been found in continuum models 
at high densities both analytically^ and numerically.— 
In the strong-coupling regime, the insulating behavior 
Nq ~ is recovered. These results are confirmed by 
GFMC (not shown), though the critical V/t is slightly 
decreased, i.e., Vc/t ^ 7. In this before with 

short-range interaction, the optimized Jastrow never be- 
haves as the potential of a CGM. 

Therefore, let us turn to the more singular interaction 
given by Eq. ([3]), leading to Vl{r) ^ — log(r). In this 
case the potential in g-space is given by ~ 1/?^ 

and we expect, similarly to what happens in the contin- 
uum for high density)^ that in the weak-coupling regime 
also Vq ~ Indeed, as shown in Fig. [THl this is the 

case for V/t < 16. Above this value, Vq becomes once 
again of the form —\og{q)/q^, just like in all previous 
examples. We note that the values of /? extracted in the 
weak coupling phase seem to be all smaller than (3* of 
the 2D Berezinskii-Kosterlitz-Thouless phase transition, 
although we can not exclude by the numerical data that 



We have shown that the long-range Jastrow wave func- 
tion gives a consistent picture of the superfluid-insulator 
transition of the bosonic Hubbard model in all spatial 
dimensions and also in the presence of long-range inter- 
action. 

In one dimension the variational results are compatible 
with a Berezinskii-Kosterlitz-Thouless phase transition 
between the quasi-long-range ordered gapless phase and 
the Mott insulator. From the point of view of the vari- 
ational wave function, the gapless phase is characterized 
by a Jastrow potential Vq ~ a/jgl, while the insulating 
one by ~ 1/g^. 

In two dimensions we have evidences of a second- 
order phase transition between a superfluid phase and 
a Mott insulator. Here, the superfluid wave function still 
has Vq ~ a/ 1 9 1, compatible with the existence of sound 
modes, while the insulating wave function is character- 
ized by a more singular Vq ~ —\og{q)/q^. This singu- 
lar behavior in the Mott phase does not change even if 
long-range interaction is considered. For instance, for a 
Coulomb interaction f2(r) ~ 1/r, the Jastrow potential in 
the superfluid phase changes into Vq ~ l/lg]^/^, compati- 
ble with the existence of 2D plasmons, yet the insulating 
wave function has still Vq ^ — \og{q)/q^ . For an interac- 
tion il(r) ~ log(r), we observe a transition between an 
algebraic long-range ordered phase^^ characterized by 
Vq ^ l/f?^, and a Mott insulating phase, once again with 
Vq log(g)/g2. 

In three dimensions, the Mott transition as revealed by 
the behavior of the Jastrow potential becomes much more 
evident. As usual, the superfluid phase has Vq ~ a/lgj. In 
this case, however, the Mott insulating wave function has 
a much more singular Vq ^ l/jql^, still with a structure 
factor that correctly behaves as Nq ^ q^ . 

This work was mainly focused on bosons, but a similar 
variational wave function can be applied also to fermionic 
models, with the same key ingredient, i.e., a long-range 
Jastrow factor that drives the metal-insulator transition, 
in spite of the uncorrelated wave function being metal- 
liCi ^^'^^ According to this picture, the Mott insulating 
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state should be produced by a sufficiently long-range Jas- 
trow potential able to bind opposite-charge fluctuations. 
Based on the bosonic results, we may argue that a Jas- 
trow potential Vq ~ l/l?!^ may work even for fermions in 
any dimensions D > 2, D — 2 playing somehow the role 
of the lowest critical dimension, Vq ^ — log(g)/g^. As 
the interaction strength is decreased, an unbinding tran- 
sition must takes place in the variational wave function, 
turning the Mott insulator into a correlated metal and 
providing a simple physical picture of the long-standing 
but still actual and very attractive phenomena which is 
the Mott transition. 

We thank important discussions with F. Alet and S. de 
Palo. This work has been partially supported by CNR- 
INFM and COFIN 2004 and 2005. 



APPENDIX A: DETAILS FOR THE 
CALCULATION OF THE STIFFNESS 

In this Appendix, we give some details for the zero- 
temperature GFMC calculation of the stiffness Dg. In 
this respect, we consider the standard Peierls substitu- 
tion and introduce an electromagnetic field in the Hamil- 
tonian by replacing the hopping term between two sites 
Ri and Rj by a suitable complex hopping)^ 



where |4'o) is the ground state of TC, with Eq energy, and 
1$) is the guiding wave function of the GFMC method. 
In the large-T limit, we have that 



ZAA] 



l ($|M/^)(vl/^|vl/o) 

T ($|*o) 



Bo) 



(A3) 



where I^^^q) and Eq are the ground state eigenfunction 
and eigenvalue of Ha, respectively. By taking the second 
derivative of Zt[A] with respect to A and then consid- 
ering A = 0, we obtain the charge stiffness (up to l/r 
corrections) . 

This quantity can be obtained by sampling sta- 
tistically the unperturbed Green's function Gx',x = 
— ^x''Hx',x/^x = Px',xbx, where Px',x is a stochastic ma- 
trix that defines the Markov chain and bx is a normal- 
ization factor. In this way, the walker \x) is distributed 
according to the variational distribution |(a;|$)p and, in 
order to obtain the true ground state, the weight G"^ must 
be considered. Then 



Zr[A] 



(A4) 



t —* t exp 




dR A{R) 



(Al) 



where the index n indicates the Markov chain iteration, 
that is defined by the transition probability Px' ,x, and 



where we can consider A{R) — {Ax,Q). Then, the first 
derivative of the energy E[A\ with respect to this field 
gives the charge current, containing both the paramag- 
netic and diamagnetic contributions. This current must 
vanish when Ax — > 0, since the Hamiltonian is real. For a 
bosonic system, the second derivative of the energy E[A\ 
represents the charge stiffness 4^ 

Within GFMC it is possible to compute the ground- 
state energy E[A\ for arbitrary time-independent static 
field Ax- Let us denote by T-La and H. the Hamiltonian in 
presence and absence of the field, respectively and con- 
sider: 



Zr[A] 



1 ($|e-"^-^|^o) 
T ($|e-^^|*o) 



(A2) 



Gl = exp 



dT'eL[x{T')] 



(A5) 



Gl[A] = G;exp<^i / dRA{R)) (A6) 



are the correcting factors of the GFMC method; eL(x) 
indicates the local energy and Rx indicates the site where 
the particle moves within the GFMC algorithm. By con- 
sidering the second derivative of Zt[A\ with respect to 
Ax, we obtain Eq. (|18p . As usual, many walkers can 
be considered with the branching technique in order to 
reduce the variance of the correcting factors G^. 
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